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This  report  is  a mathematical  study  of  time-domain  Wiener  adaptive 
multichannel  filtering  using  distributed  signal  models.  It  presents  derivations 
for  two  Wiener  adaptive  algorithms  and  describes  the  basic  technique  for  im- 
plementing a Wiener  adaptive  processor  with  directionally  distributed  signal 
models. 


In  order  to  obtain  the  signal-model  crosscorrelation  functions  needed 
in  die  Wiener  adaptive  algorithms,  the  probability  density  function  p(r)  for 
die  time  delay  between  a reference  sensor  and  the.  individual  channels  must 
be  specified.  Analytic  derivations  of  this  function  are  presented  for  inverse 
velocity  space  models,  distributed  ring  models,  and  velocity-azimuth  space 
models.  These  derivations  arc  of  interest  in  their  own  right:  they  are  useful 

in  specifying  two-channel  crosscor relation  functions  for  various  directional 
energy  distributions. 
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SECTION  I 


INTRODUCTION 

This  report  presents  the  results  of  a mathematic  study  directed 
toward  implementation  of  Wiener  adaptive  multichannel  filtering  with  dis- 
tributed signal  models. 

Section  II  discusses  the  reasons  for  using  this  form  of  Wiener 
filtering,  derives  both  an  unconstrained  Wiener  adaptive  algorithm  and  a 
Wiener  adaptive  algorithm  subject  to  the  same  unity-response  constraints 
as  specified  in  the  maximum  likelihood  adaptive  algorithm,  and  describes 
a technique  for  modeling  the  signal  crosscorrelation  functions  needed  in 
the  Wiener  adaptive  algorithms. 

Sections  III  through  V give  analytic  derivations  for  the  time-lag 
probability  density  function  corresponding  to  specific  distributed  signal 
models.  The  time-lag  probability  density  function  p(r)  for  the  time  delay 
between  the  reference  sensor  and  the  individual  channels  is  necessary  in 
order  to  obtain  the  signal-model  crosscorrelation  functions  in  the  Wiener 
adaptive  algorithms.  Section  III  features  inverse  velocity  space  models, 
Section  IV  distributed  ring  models,  and  Section  V velocity -azimuth  space 
models. 

Section  VI  summarizes  the  results  of  diis  mathematical  study. 
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SECTION  II 


IMPLEMENTATION  OF  WIENER  ADAPTIVE  FILTERING 
WITH  DISTRIBUTED  SIGNAL  MODELS 

A.  DISCUSSION 

A previous  adaptive -filter ing  study  (Barnard  and  O'Brien,  1974)  em- 
ployed a maximum  likelihood  adaptive  algorithm  where  the  adaptive  filter  out- 
put y(t)  was  formed  using  the  equation 

y(t)  = xta 

and  the  adaptive  filter  vector  A was  updated  according  to  the  equation 

A(t  +At)  - A(t)  = Zn(t)*1  A(X-X)  = 2/l(t)y(t)(X-X)  , 

where  the  filter  weight  vector  A , the  data  vector  X , and  the  beamsteer  out 
put  vector  X are,  respectively, 
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jU(t)  is  a scalar  quantity  controlling  the  convergence  rate.  The  superscript  T 
denotes  vector  transpose.  M is  the  number  of  channels,  and  2N+1  is  the 
total  length  of  the  adaptive  filter  in  points.  After  time-shifting  the  input  chan- 
nels to  align  energy  from  the  specified  look  direction,  the  quantities  x.(t-j)  in- 
side the  data  vector  X are  the  values  for  channel  i at  time  t-jAt  , where 
At  is  the  sampling  interval.  Inside  each  subvector  of  the  beamsteer  output 
vector  X , all  vector  components  are  identical  and  equal  to  the  beamsteer  out- 
put 

M 

- TT  X) 

i=  1 


at  time  t - jAt. 


When  operating  against  background  noise  from  the  Alaska  Long  Period 
Array  (ALP A),  the  maximum  likelihood  adaptive  algorithm  (as  implemented  in 
the  study  mentioned)  functioned  well  and  produced  6 dB  signal-to-noise-ratio 
imorovements  with  suitably  chosen  operating  parameters.  In  interfering-event 
simulations  using  ALPA  data,  the  same  algorithm  produced  significant  detec- 
tion gains  but  distorted  the  on-azimuth  signal  waveform  considerably.  The 
distortion  was  principally  a result  of  the  maximum  likelihood  design  goal 
(minimization  of  the  filter  output  power  y^t)  subject  to  unity-response  con- 
straints in  the  beam  look  direction).  In  conformity  with  this  design  goal,  the 
adaptive  processor  attempted  to  create  from  the  data  sample  containing  the  off- 
azimuth  interfering  event  a filter  output  180°  out  of  phase  with  the  data  sample 
containing  the  on-azimuth  event  and  was  partially  successful  in  producing 
mutual  cancellation  of  the  two  events.  A possible  remedy  for  this  problem  is 
a Wiener  adaptive  filtering  algorithm,  where  the  design  goal  is  to  minimize 
the  mean  square  difference  between  the  adaptive  filter  output  and  the  on-azimuth 
signal.  Conceivably,  such  an  algorithm  might  even  provide  substantial  detec- 
tion improvement  over  the  maximum  likelihood  algorithm  actually  employed. 
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The  subject  of  this  section  is  the  basic  technique  for  implementing  a 
Wiener  adaptive  processor  using  distributed  signal  models.  Discussion  of 
specific  distributed  signal  models  is  confined  to  later  sections  of  this  report. 
Subsection  B in  this  section  presents  derivations  for  both  an  unconstrained 


Wiener  adaptive  algorithm  and  a Wiener  adaptive  algorithm  subject  to  the  same 
unity-response  constraints  in  the  look  direction  as  specified  in  the  maximum 
likelihood  adaptive  algorithm.  Both  of  the  derived  algorithms  require  some 
method  for  estimating  the  c ros  sco r relation  functions  between  the  signal  and 
the  charnels  entering  the  adaptive  beamformer.  Subsection  C describes  a 
technique  for  modeling  the  crosscorrelation  functions  needed  in  the  Wiener 


adaptive  algorithm. 


B.  DEFLATION  OF  WIENER  ADAPTIVE  FILTERING  ALGORITHMS 

In  the  case  of  the  unconstrained  Wiener  adaptive  filtering  algorithm, 

the  derivation  is  relatively  straightforward.  The  error  e(t)  in  estimating 

T 

the  signal  s(t)  is  € (t)  = s(t)  - A X and  the  mean  square  error  is 
€Z(t)  = [s ( t ) - ATx][s(t)  - X1  a] 

= sZ(t)  - 2A7  X s ( t ) + A1  XX1  A . 


When  the  Widrow  approximation  XX  ~ XX  is  used  in  this  expression,  the 
gradient  of  the  mean  square  error  is 


2XXTA  - 2V  , 
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corresponding  to  a white  frequency  response  in  the  look  direction.  Symbolically, 
the  constraint  matrix  equation  may  be  written  CA=D,  where  C is  the  2N  + 1 by 
M(2Nt  1)  constraint  matrix  and  A is  the  M\2N+  1) -dimensional  adaptive  filter 
vector.  The  adaptive  filter  must  move  in  the  direction  opposite  the  gradient 


Ve  (t)  - A (D  - CA) 


} 


T 

XX  A 


2 V + 


T 

C A 


where  A is  the  Lagrangian  multiplier  vector 


A = 


X ( ~N) 

X(-l) 

A(0) 

X(l) 


LX(N)  J . 


The  update  equation  is 


A(t+  At) 
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/'(t) 
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and  the  Lagrangian  multiplier  vector  is  determined  by  premultiplying  this 
equation  by  the  constraint  matrix  C: 


D = CA  (t  + At)  = CA  - ji(t)  |2CXXTA  - 2CV  + CCTa]  . 


Rearrangement  of  this  equation  yields 

jl(t)  CCT\  = [cA-dJ  + 2/Ll (t)  jcv  - CXXTaJ 


or 

A=  2^CCT^  [cV-CXX^a], 
so  that  the  adaptive  filter  update  equation  is 


A(t+  At)  - A(t)  = 2/1  (t)  [CT  (CCT)  C - l][xXTA  - v]  . 


The  matrix  GC  is  the  (2N+1)  by  (2N+1)  diagonal  matrix 


M 0 


0 MO. 


'•  0 M 0 


1)  M 


and  the  matrix  C 


T 


(cCT)  C i 


s the  M(2N+1)  by  M(2N+1)  matrix 


. i 
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where  each  of  the  2Nfl  submatrices  is  an  M by  M matrix  with  each 

element  equal  to  one.  Premultiplication  of  any  M (2N+1  )-dimensional  column 
vector  by  the  matrix  C produces  a column  vector  where  all  M elements  in 
each  of  the  2N+1  subvectors  are  identical  and  equal  to  the  average  of  the  corres- 
ponding elements  in  the  original  column  vector.  As  a result,  the  adaptive  filter 
update  equation  reduces  to 

A(t  + At)  - A(t)  - 2M(t)  Jy(t)  (X  - X)  - (V  - V)J  , 
where  each  element  of  the  vector 
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is  equal  to 

M 

^(t_j)  = h.  X)  xi(t-j)  s(t) 

i=  1 

If  the  adaptive  filter  vector  A initially  satisfies  the  constraint  matrix 
equation  CA  = D,  the  adaptive  update  equation  produces  adaptive  filter  vectors 
A which  also  satisfy  the  constraint  equation.  Because  of  computer  roundoff 
error,  however,  the  filter  vector  A must  periodically  be  adjusted  to  satisfy 
the  constraints. 

C.  ESTIMATION  OF  THE  CROSSCORRELATION  VECTOR 

Since  seismic  adaptive  processing  must  be  performed  in  the  time  do- 
main to  obtain  continuous  time  traces,  the  technique  described  in  this  sub- 
section computes  the  signal  model  crosscorrelation  function  components 
x.(t.-j)  s ( t)  of  the  vector  V = Xs(t)  directly  in  the  time  domain  as  a function 
of  the  time  lag  r between  the  signal  as  seen  by  the  reference  sensor  and  the 
signal  as  seen  hy  the  i-th  sensor. 

If  the  power  spectrum  of  the  signal  were  white,  the  required  cross- 
correlation function  would  simply  be  the  probability  density  function  p(r) 
for  the  signal's  time  lag  t between  the  reference  sensor  and  the  i-th  sensor. 
For  a given  distance  d between  the  two  sensors,  the  probability  density  func- 
tion for  T = d cos  0/  V depends  on  the  joint  probability  distribution  for  0 and 
V , where  0 is  the  arrival  angle  of  the  signal  and  V is  the  apparent  velocity 
of  propagation  (see  Figure  II - 1 ) . 

In  order  to  reflect  the  signal  power  spectrum,  the  crosscorrelation 
function  components  x.(t-j)  s ( t)  are  obtained  by  convolving  the  probability 
density  function  p(r)  for  the  specified  distributed  signal  model  with  the  auto- 
correlation function  corresponding  to  the  estimated  signal  power  spectrum. 
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is,  for  incoming  energy,  the  arrival  angle  (in  the  plane  of  the 
array)  relative  to  the  line  joining  the  two  sensors 

is  the  apparent  velocity  of  propagation  with  respect  to  the  plane 
of  the  array 

is  the  distance  between  the  two  sensors 

is  the  time  lag  of  incoming  energy  between  the  two  sensors 

d cos  0 

T = — 


FIGURE  II- 1 

DIAGRAMMATIC  REPRESENTATION  OF  THE  FACTORS  DETERMINING 
THE  TIME  LAG  OF  THE  SIGNAL  BETWEEN  TWO  SENSORS 


II- 9 


By  choosing  the  power  spectrum  of  the  data  as  the  model  for  the  signal  spectrum, 
the  adaptive  b amformer  can  be  made  to  function  solely  as  a directional-filtering 
processor. 

For  array  sensors  with  well-equalized  signal  outputs,  the  average  auto- 
correlation function  across  channels  can  serve  as  a model  for  the  signal  auto- 
correlation function: 

M M 

s(t-k)  s(t)  = ~ ^ x (t-k)  x (t)=ri  ^ ' x (t-k)  x (t) 

m=l  m=l 

If  the  sensor  responses  differ  by  a simpl;  scale  factor,  then  an  amplitude 
correction  can  be  made1 


M 


5.  (t-k)  s . ( t)  = 1-  ^ x (t-k)  x (t) 
1 t M m m 


m=  1 


where  (0)  and  (P  (0)  denote  the  average  power  for  the  channels  denoted 
ii  mm 

by  the  subscripts  i and  m,  respectively.  (Since  the  crosscorrelation  function 

x.(t-j)  s(t)  is  being  estimated  and  the  strength  of  the  i-th  sensor  affects  only 

the  t^rni  x.(t-j),  only  an  amplitude  correction  is  made  instead  of  a power 

correction.  ) For  good  results  in  eliminating  the  variations  in  sensor  response, 

the  power  ratios  (/?..( 0 )/  a)  (0)  need  to  be  reasonably  stable. 
r n mm 

In  the  later  sections  of  this  report,  formulas  for  the  probability  density 
function  p (r)  will  provide  expressions  for  p(r)  as  a continuous  function  of 
the  time  lag  T . The  n -nzerc  values  of  p(r)  lie  within  the  interval 


d 

max 


V . 
min 


d 

<T  <_  maX 

” “ V . 

mm 


> 
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where  d is  the  maximum  distance  to  the  reference  sensor  and  V^ 

is  the  minimum  apparent  propagation  velocity  across  the  plane  of  the  array, 

so  that  the  frequency  bandwidth  for  p(r)  is  infinite.  On  the  other  hand,  the 
autocorrelation  function  to  be  convolved  with  p(r)  in  estimating  *.(t-j)  s(t) 
is  available  only  at  integer  multiples  of  the  sampling  interval  At.  In  well- 
designed  processing  systems,  spectra  of  signals  and  noise  lie  nredominantly 

within  the  frequency  band 

_ W_  < f < w 

where  W = 1/At  is  the  bandwidth  corresponding  to  the  sampling  interval 
At,  and  it  is  not  a bad  approximation  to  assume  that  the  spectrum  of  the  data 
lies  entirely  within  the  Land  indicated.  If  the  unmodified  probability  density 
function  p(T)  were  sampled  and  then  convolved  with  the  available  values  ot 
the  estimated  signal  autocorrelation  function,  then  the  frequency  components 
of  p(T)  outside  the  sampling  bandwidth  would,  in  effect,  alias  back  into  that 
bandwidth.  A remedy  for  this  problem  is  to  compute  a continuous  convolution 

of  p (t)  with 


sine  Wr  = 


sin  TT  WT 


ttwt 


if  7=0 


it  r =j£  0 


at  integer  multiplies  of  the  sampling  interval  At  by  numerical  integration.  The 
number  of  computed  values  should  be  sufficient  to  extend  at  least  beyond  all  lags 
T for  which  p(T)  is  nonzero.  Convolution  with  sine  tv 7 removes  from  p(t) 
all  spectral  components  outside  the  sampling  bandwidth.  If  an  inlinitc -length 
version  of  this  bandlimited  and  sampled  version  of  p(T)  were  convolved  with 
an  infinite-length  sampled  autocorrelation  function  estimate,  the  result  would 
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be  the  true  sampled  convolution  of  the  probability  density  function  p(r)  with 
the  bandlimited  autocorrelation  function  estimate.  In  practice,  of  course, 
both  of  these  functions  must  be  truncated  to  a finite  length.  If  the  point  length 
is  2L,+  1 for  the  modified  version  of  p(T)  and  2(N+Li)  +1  for  the  autocorrelation 
function  estimate,  then  the  Fourier  transform  of  the  (2 N+l )-point-long  convolu- 
tion output  is  the  desired  signal-model  crosspower  spectrum  convolved  with 
Jsin  IT  (2  L+l ) f A t J / f^sin  77  f A tj  , which  is  the  Fourier  transform  of  the 
(2  L+l  )-point-long  sampling  function  for  the  modified  probability  density  function 
p(r).  If  the  number  of  points  2L+  1 is  sufficiently  large  and  the  estimated  signal- 
model  power  spectrum  is  sufficiently  smooth,  the  signal-model  ci  osscorrelation 

function  estimates  for  x (t-j)  s(t)  will  be  relatively  free  of  spectral  window 

i 

effects. 

In  the  following  sections,  analytic  derivations  are  given  for  the  probabil- 
ity function  p(r)  corresponding  to  specific  distributed  signal  models.  These 
derivations  generally  are  obtained  using  a transformation  from  the  specified 
velocity-azimuth  probability  distribution  to  a new  two-dimensional  space  where 
the  time  lag  r is  one  of  the  coordinat'  s. 


. i 
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SECTION  III 

INVERSE  VELOCITY  SPACE  MODELS 


A.  BASIC  APPROACH 

In  this  section,  the  signal-model  probability  distributions  are 
specified  as  a function  of  the  two-dimensional  inverse  velocity  vector 

— ► 

— V 

u = — » 

— * — > 

v • v 


where  v is  the  incoming  energy's  apparent  velocity  in  the  plane  of  the  array. 


The  coordinate  system  axes  are  oriented  so  that  the  ux  coordinate  points  in  the 


same  direction  as  the  line  from  the  i-th  sensor  to  the  reference  sensor 
(as  shown  in  Figure  III- 1).  The  time  lag  t corresponding  to  a particular  value 
, — * 

of  U IS 


T = d u , 

X 

where  d is  the  distance  between  the  two  sensors,  and 


~dT  ~ d 

The  probability  density  function  p (tq)  is  determined  by  taking  the  limit  as 
AT  approaches  zero  for  the  quantity  (1/  At)  times  the  probability  between  T 

and  t + At: 


J<k 


is  the  vector  location  of  the  i-th  sensor 

is  the  vector  location  of  the  reference  sensor 
r 

d is  the  distance,  between  the  two  sensors 


* V 

u = — 

v"-  V 

u is  the  inverse  velocity  vector  corresponding  to 
the  apparent  velocity  across  the  array 

V is  the  incoming  energy's  apparent  velocity  in 
the  plane  of  the  array 


FIGURE  m-1 

THE  BASIC  INVERSE  VELOCITY  SPACE  MODEL 


lim  | P<  r0  * ro  + Ar)  P K’o-VAV  AUx] 
P{To}  = At-VO)  Tt  = At 


1 p (u  ) , u du  A a 

lim  ImL/L  xo  Vi  vJ  x 

Ar-*0  AT 


00 

j p[("x»c’  Uv]dU: 


B.  DISK  MODEL 


For  this  model,  the  probability  density  is  uniform  on  the  disk  [u|  < 
where  u = 1/V  and  V is  the  minimum  apparent  velocity  across  the 


array: 


-L  « v 

2 ’ max 

i r u 

max 


p(u  , u ) 

^ y 


0 if 


M >u 


max. 


The  probability  density  function  p(r)  for  the  time  lag  t is 


2„  / u 


77u  d 

max 


77  u d 
max 


for  |r|  < d/V  aaci  zero  elsewhere.  Figure  111-2  displays  p(-r)  for  this  model, 
whose  principal  application  arises  in  connection  with  infinite -veloc ity  P waves. 


C.  ANNULAR  MODEL 

In  this  model,  the  probability  density  is  uniform  on  the  annulus 

u . <|uf<  u , where  u = 1 / V .,u  . = 1/V  , and  V . and  V 

min  max  max  min  min  max  min  max 

are  the  minimum  and  maximum  apparent  velocities,  respectively,  across  the 
array: 


p(u  , u ) 
x y 


if 


u 

min 


< u< 


if  fu|  < u . or 
min 


u 

max 


> u 

max 


The  probability  density  function  p(r)  for  the  time  lag  r is 


if  u < u 

x ~ min 


if  u < u <u 

min  — x~  max 


and  zero  for  |r|>d/V  . Figure  III  - 3 presents  p(r)  for  the  annular  model 

at  nine  values  of  V . / V 

min  max 

D.  SECTOR-OF -ANNULUS  MODEL 

In  this  model,  the  probability  density  is  uniform  on  the  sector  of  the 

annulus  u . < u|  < u between  0 . and  n , which  are  the  mini- 

mm  max  min  max 

mum  and  maximum  arrival  angles  relative  to  the  line  joining  the  i-th  sensor 
and  the  reference  sensor.  That  is  to  say, 


^ 2 

1(0  - 0 . )(u2  - u2  . ) 

I max  min  max  min 


if  u . < 

min 


u 

max 


and  f) 

min 


< 0 < Q 


max 


otherwise. 


Breaking  up  the  angular  range  0 . < 0 < 0 into  sections  bounded  by  the 

min  max  ’ 

transition  angles  encountered  as  0 sweeps  from  0 . to  0 con- 

min  max 

siderably  simplifies  the  determination  of  p (r).  If  0 . and  Q are  speci- 

min  max 

fied  so  that  -2  7r  £_0  . < 2 n,  -2tt  < 0 <^2  tt  , and  Q . < 0 < <9  . d-  2 tt , 

min  max  min  max  ■*-  "min 

then  there  are  no  more  than  three  sections,  each  lying  either  in  the  first  and 
second  quadrants  or  in  the  third  and  fourth  quadrants.  Sections  can  be  mapped 
into  the  first  and  second  quadrants  in  such  a way  that  the  probability  density  as 
a function  of  r is  unaffected.  Figure  III-4  depicts  a procedure  for  determining 


the  minimum  angle  a and  the  maximum  angle /f  of  each  section  after  the  appro- 
priate mapping  onto  the  first  and  second  quadrants.  Once  the  subdivision  of  the 

angular  range  0 < a <-  a and  the  mapping  onto  the  first  and  second  quad- 

min  u "max 

rants  have  been  completed,  the  probability  density  function  p(r)  is  determined  by 
adding  the  contributions  from  the  appropriate  sections  at  the  time  lag  T.  If  the 
minimum  angle  Of  and  the  maximum  angle  P after  mapping  onto  the  first  and 
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FIGURE  III- 3 

TIME-LAG  PROBABILITY  DENSITY  FUNCTION 
FOR  THE  INVERSE  VELOCITY  SPACE  ANNULAR  MODEL 
AT  NINE  VALUES  OF  V^/V^^ 
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N = 0 

M = 1 


FIGURE  III -4 

PROCEDURE  FOR  MAPPING  ANGLES  BETWEEN 

i AND  a ONTO  FIRST  AND  SECOND  QUADRANTS 
min  max 


second  quadrants  both  lie  in  the  same  quadrant  and  the  annulus  between  u 

min 

and  u is  thick  enough,  there  will  be  a ranee  of  u for  which  the  extremal 
max  x 

values  of  u are  located  on  the  section  boundaries  u = (u  cos  /?  , u sin/?)  and 

u = (u  cos  (V  , u sin  a ).  On  the  other  hand,  if  a lies  in  the  first  quadrant  and  /? 

in  the  second  quadrant  or  if  the  annulus  is  thin  enough,  there  will  be  a range  of 

u for  which  the  extremal  values  of  u lie  on  the  circles  Itfl  = u and  |u*|=  u 

x y max  min 

Figure  III  - 5 illustrates  these  two  possibilities  when  cr  and  [3  both  lie  in  the  first 

quadrant.  The  case  shown  in  Figure  III- 5 ( a.)  corresponds  to  the  condition  u cos R 

max  r 

^u  . cos  a,  while  the  case  shown  in  Figure  III- 5(b)  corresponds  to  the  condition 
min  ° r 

u . cos  a < u cos  B . 
min  max 

On  the  line  u = (u  cos  y,  u sin  y),  u - t / d and  u = r tan  7 / d,  so 

x y 

that  u =rtan  a/ d or  u = t tan  Bid  whenever  a section  boundary  corresponds 

y y 

to  the  extremal  value  for  u^  at  the  time  lag  T . For  any  given  pair  of  angles 
n and/3  , there  are  three  possible  u^  intervals  where  the  formula  for  p(r)  is 
different.  In  the  first  and  last  of  these  intervals,  one  extremal  value  of  u is 

y 

on  a section  boundary,  while  the  other  is  on  a circular  arc.  In  the  middle  inter- 
val, both  extremal  values  lie  either  on  circular  arcs  or  on  section  boundaries. 

Let 

2 


K = 


d(  0 - 0 . )(u2  - u2  . ) 

max  min  ' max  min' 


/ 1 

1 \ 

®max"  ^\r\ir)  ( y2  . 

\ mm 

v2  ) 

max/ 

If  ff  /2  < B , then  the  first  interval  is  u cos R < u < min(u  cos  B , 

max  ' ~ v 

cos  o ) and 


mm 


u 


max 
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(a)  No  Time-  .Lag  Overlap  Between  Section  Boundaries 


(b)  Time-Lag  Overlap  Between  Section  Boundaries 

FIGURE  III -5 

THE  TWO  FIRST-QUADRANT  POSSIBILITIES  AFFECTING  THE  FORMULAS 
FOR  THE  TIME-LAG  PROBABILITY  DENSITY  FUNCTION 
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P(  T ) 


7 2 

/u 

2 

- u 

Ttan/3 

>/  max 

X 
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|_  mm  V 


for  (d  cos  (3/V  . ) < t <_ min  [(d  cos  /3/V  ),  (d  cos  « / V )].  If/3  < v /2 


mm’  max 

then  the  first  interval  is  u cos  (3  <^n  < min  (u  cos  [3  , u cos  n ) and 

i « v m ax  ^ * ' ' 


min 


x 


min- 
min 


P(  r ) = K 


= K 


’ tan/3 


-/ 


2 2 

u . - u 

min  x 


rtan/3 

__ 


V 


max 


f 


for  (d  cos  /3/V  ) <_  r < min  [(dcos/3/V  ),  (dcosar/V  )].  If 

' * *-^1  n v min  iTlrlX 


rnax 


max 


/3  - n /2,  the  first  interval  is  of  zero  length,  and  the  formula  for  the  second 

interval  should  be  used  at  T - 0.  If  max  (u  cos  (3  , u cos  (3  ) <_u  <. 

max  min  x 

min  (u  cos  a,  u . cos  a),  then  the  extremal  values  of  u in  the  second 
max  min  y 

interval  are  on  circular  arcs  and 


P(r)  = 


r / 2 

2 

nr 

2 ~\ 

/ u 

- u 

- 

Vu 

u 

max 

X 

v min 

X 

/ V . T\Z 

I mm  I 

\ “ d / 


V 


max 


for  max|(d  cos/3/V  . ),  (d  cos  /3/V  )]  _£.  r < min  [ (d  cos  a /V  . ), 

1 mm  max  mm 

(d  cos  «/V  )1  • If  min  (u  cos  a , u cos  a ) <_u  < max  (u  cos  (3, 

max  max  mm  x—  max 

U cos  B),  then  the  extremal  values  of  u in  the  second  interval  are  on  section 
min  ^ y 

boundaries  and 

p(T)  = — J T — |tan  (3  - tan  o;  I 
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■fc. 


for  min  I (cl  cos«/V  . ),  (d  toscv/V  )]  < r < max  [ (cl  cos  (3  / V ), 

min  max  mm 

d cos/?  / V )1.  If  the  limits  for  the  second  interval  are  equal,  then  only 

max 

the  formulas  for  the  first  and  third  intervals  should  be  used.  If  7t/2<cy  , 

then  the  third  interval  is  max  (u  . cos  U , u cos  cv)  u u cos  « and 

min  max  x nun 


r ) = K — ~ 


rtan« 

~ d 


ru  . - u 

min  > 


= K 


T tanev 


'V  r\ 
max 


for  max  I (d  cos  /?  / V ),  (d  cos  cv  / V . )1  It  <(d  cos  (»  / V ).  If  a <7r/  2, 

1 max  min  max 

then  the  third  interval  is  max  (u  cos  13  , u . cos  « ) < u < u cos  « and 

max  min  *"  x “ max 

% 1/2  2 rtann;  "| 

p(  r)  K /u  - u - t 

[_v  max  x d 

r M/  r \" "2 

T r 1 /,  I min  1 r lan  (v 

K — V1  • V~T -7  - — d- 

min  J 

for  max  [ (d  cos  a/  V ),  (d  cos  a/V  )]  ^ t<  (d  cos  « /V  . ).  If  « = tt  / 2, 

then  the  third  interval  is  of  zero  length,  and  the  formula  for  the  second  interval 

should  be  used  at  r=  0. 

Figures  III  - 6 and  III  - 7 portray  the  time-lag  probability  density  function 

p(r)  at  V / V =0.1  and  V . /v  = 0.9,  respectively,  for  vari- 

min  max  nun  max 

ous  30  sectors. 
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TIME-LAG  PROBABILITY  DENSITY  FUNCTION  FOR  TIIF  INVERSE 
VELOCITY  SPACE  SECTOR  -OF-ANNULUS  MODEL 

AT  V /V  - 0.  1 FOR  VARIOUS  30°  SECTORS 
min 


max 


SECTION  IV 

DISTRIBUTED  RING  MODEL 


In  this  form  of  distributed  signal  model,  all  incoming  energy  is  con- 
centrated at  a single  apparent  velocity  V with  respect  to  the  plane  of  the 
array.  The  incoming  energy,  however,  is  distributed  over  a range  of  arrival 
angles  with  respect  to  the  line  joining  the  i-th  sensor  to  the  reference  sensor. 

If  the  probability  density  function  p q(0)  is  known  between  the  arrival  angles 
. tt  and  n,  the  probability  density  function  P(t)  for  the  time  lag  T can  be  easily 
determined.  First,  a transformation  from  velocity-azimuth  space  to  delay- 
azimuth  space  is  made  with  the  formula 

_ d cos  9 , 

V 

where  d is  the  distance  between  the  two  sensors.  The  line  V=  constant  is 
transformed  into  a cosine  curve  (as  shown  in  Figure  IV- 1).  For  a given  time 
lag  r , the  corresponding  angles  are  0 - cos  (V  r / d)  and  0 - - cos  (V  T /d). 


When  any  time  lag  Tq  is  less  than  -d/V,  the  probability  P(t  £_Tq)  that 
T is  not  greater  than  T is  the  constant  value  zero.  Correspondingly,  when 
any  time  lag  T is  greater  than  d/V,  the  probability  P(t  Tq)  is  one.  Thus 
the  probability  density  function  p(t)  is  zero  for  |r|>d/V.  Between 
T - - d/V  and  T = d/V, 


P(TiT0) 


( 

-i  (Vr«A 

% 

cos  y— ; 

" Fd 

cos  \ d / 
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- 

FIGURE  IV-1 

TRANSFORMATION  FROM  VE  LOCITY- AZIMUTH  ?°ACE 
TO  DELAY- AZIMUTH  SPACE  FOR  DISTRIBUTED  RING  MODEL 


where  F q{0)  is  the  cumulative  distribution  function  for  the  ariival  angle  6 . 
Thus,  for  | r I < d/V,  the  time-lag  probability  density  function  is 


dP(  T < T ) 

o 


If  the  arrival  angle  probability  density  function  p f)(0)  is  nonzero  when  0 is  ti 
integer  multiple  of  TT , there  will  be  a discontinuity  in  p(r)  at  M = d/V:  a s | r | 
approaches  d/V  from  smaller  absolute  values,  p(r)  goes  to  infinity;  for  all 
values  of  I r | greater  than  d/V,  on  the  other  hand,  p(t)  is  zero. 


Figure  IV -2  depicts,  for  select  d 30  sectors,  the  probability  density 
function  p(r)  corresponding  to  the  uniform  angular  probability  density  func- 


p (9)  - 

V 


if  0 . < 0 <9 

0 -a  min  max 

max  min 


otherwise. 


In  the  case  where  0 = -10°  and  ()  = 20°,  the  probability  density 

min  max 

approaches  oo  as  T approaches  d/V  from  the  left. 

Figure  IV- 1 pictures  p(t)  when  the  incoming  energy  is  equally  likely 
at  all  possible  arrival  azimuths.  The  time-lag  probability  function  approaches 
oo  as  T approaches  ± d/ V moving  outward  from  T-  0.  This  figure  illustrates 
the  limiting  case  for  the  inverse  velocity  space  model  when  Vmax 

approaches  one  (compare  with  Figures  III  - 2 and  III-3). 

Figure  IV-4  represents  the  probability  density  function  p(T)  for  the 
time  hag  T when  the  arrival-angle  probability  density  is  normally  distributed 
about  some  angle  0 : 


(o-o0) 


po{0)  = ~r — ° 

For  this  figure,  the  standard  deviation  CT  is  15  . Strictly  speaking,  there  is 
an  infinite  discontinuity  at  T = ± d/V  in  all  cases,  but  the  probability  associated 
with  these  discontinuities  is  sometimes  minuscule  and  the  discontinuity  for  such 


cases  is  not  shown  in  the  figure. 


IV  - 4 


FIGURE  IV - 2 

TIME-LAG  PROBABILITY  DENSITY  FUNCTION  FOR  THE  DISTRIBUTED  RING 
MODEL  WITH  A UNIFORM  DISTRIBUTION  OVER  VARIOUS  30°  SECTORS 


FIGURE  IV -4 

TIME-LAG  PROBABILITY  DENSITY  FUNCTION ^ FOR ^ 
MODEL  WITH  A NORMAL  DISTRIBUTION  (a  - 1 5 ) PEAKING  A J 

VARIOUS  ANGLES 


SECTION  V 

VELOCITY-AZIMUTH  SPACE  MODELS 


A. 


B/  SIC  APPROACH 
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In  this  section,  the  signal-model  probability  distributions  are  specified 
as  a function  of  the  apparent  velocity  V in  the  plane  of  the  array  and  tin* 
arrival  angle  0 with  respect  to  the  line  joining  the  i-th  sensor  to  the  refer- 
ence sensor.  The  cumulative  distribution  functh  n F(  r ) for  the?  time  lag 

T , i.  e.  , the  probability  P(r<  r ) that  r is  less  than  or  equal  to  some  arbi- 
o o 

trary  time  lag  r , is  obtained  by  integrating  over  all  points  (V,  6 ) corres- 
ponding to  time  lags  less  than  or  equal  to  r : 


F(r  ) 
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P(r_<  r ) - 


[J 

T^o 


p(V,  0)  d 0 dv, 


where  p(V,0)  is  the  probability  density  at  the  point  (V,fl).  The  transformation 
of  variables 


T = 


d cos  8 


0 = 0 


1 


T 

■ 't 

i 
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facilitates  the  evaluation  of  the  double  integral  above.  If  the  probability  density 

is  nonzero  only  inside  the  rectangle  specified  by  V . < V<  V and 

7 min  max 

- it < 7r,Figur c V-l  illustrates  the  corresponding  area  in  ( r,  0)  space.  The 
shaded  area  corresponds  to  the  specified  rectangle  in  (V,0)  space.  After 
transformation,  the  cumulative  distribution  function  F(  r ) is 
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TR  ANSFORMATI' 


■vhei  '•  J{t,0)  is  the  Jacobian 
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The  probability  density  function  for  the  time  lag  r is 


P (t) 


dF(r) 
cl  r 
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P (v,0) 


For  bivariate  polynomial  probability  densities 
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p(v,»)  = 22  22 


a..  vV1 


ij 
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i=0  i=0 


cos  ()  d 0 


the  integrand  is  of  the  form 


M N 

^ ^ ^ ^ b_  (cos 0)  1+10^ 

i=0  j=  0 


and  can  be  evaluated  analytically. 


B.  UNIFORM  RECTANGULAR  DISTRIBUTION 

In  this  model,  the  probability  density  is  uniform  over  a rectangle 
in  velocity-azimuth  space: 


^Vmax'Vmin^  0max  ^min) 


p;V,  0)  = 
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if  V . <V<  V 

min  max 

and  0 ■ <0<O 

inin^"  "nax 


otherwi  se. 


The  time-Ug  probability  density  function  is 


“ |t(t  ( V - V . )(  o - 0~J 

max  mm  max  min 
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cos  0 d 0 
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d sin  6 


^T^Vmax  Vmin^  flmax  ^min^ 


0(t) 


If  the  shaded  areas  above  #max  and  below  0^^  in  Figure  V-l  are  removed, 
then  the  rema  ning  shaded  area  determines  the  limits  of  integration  for  0 at 
any  particular  time  lag  r . Since  the  angles  0 along  the  curved  boundaries 
in  Figure  V-l  are  of  the  form  cos  (Vr/d)  and  -cos  (Vr/d),  the  following 
identities  are  useful  in  the  integral  evaluation: 
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As  in  the  case  of  the  Section  III  sector-of-annulus  model,  subdividing  the 

angular  range  0 . <0<0  into  sections  bounded  by  the  multiple -of- tt 

min  max 

transition  angles  encountered  as  0 progresses  from  0^  to  0^  facilitates 
the  determination  of  p(r).  As  before,  sections  can  be  mapped  into  the  first 
and  second  quadrants  so  that  the  probability  density  as  a function  of  r is  un- 
affected. The  same  procedure  outlined  in  Figure  III-4  is  also  sufficient  to 
ascertain  the  minimum  angle  « and  the  maximum  angle  /3  of  each  section 
after  mapping  onto  the  first  and  second  quadrants.  Addition  of  the  contributions 
from  the  appropriate  sections  at  the  time  lag  r again  determines  the  probabil- 
ity density  function  p(  r).  If  0^  and  0^  are  specified  so  that  - 2 tt  < »min<  2 tt, 

<2tt,  and  0 . <0  < 0 . + 2tt,  then  there  are  no  more  than 

c < “max — min  max  mm 

three  such  sections. 

The  extremal  values  of  0 for  any  particular  time  lag  r lie  either  on 
horizontal  section  boundaries  or  on  arc-cosine  curves.  The  boundaries  in 
the  rectangular  velocity-azimuth  model  correspond  precisely  to  the  boundaries 
of  the  inverse-velocity  sector-of-annulus  model  in  Section  III.  The  difference 
between  the  two  models  is  the  probability  density  within  the  boundaries:  if  the 
uniform  rectangular  velocity-azimuth  model  were  expressed  in  terms  of  the 
inverse  velocity  model,  the  probability  density  would  rise  as  the  vector  magni- 
tude [u|  decreases.  As  in  the  case  of  the  sector-of-annulus  model,  there  are, 
for  any  give  pair  of  angles  « and  (3  , three  possible  r intervals  where  the  formula 
for  p(r)  is  different.  In  the  first  and  last  of  these  intervals,  one  extremal  value 
of  0 is  on  a horizontal  section  boundary,  while  the  other  is  on  an  arc-cosine 
curve.  In  the  middle  interval,  both  extremal  values  lie  either  on  arc -cosine 
curves  or  on  horizontal  section  boundaries. 


V-5 


If  7T / 2 < B . then  the  first  interval  is  (cl  cos  /3  / V . ) < r< 
^ min  — “ 

min  I (d  cos  B/V  ),  (d  cos  cv/V  . )1  and 
1 max  mm  J 


d | sinfi  - sin  cos 


V . T 
min 


r r (V  - V . )(  0 - 0 . ) 

max  min  max  min 

/a  - (V  . r )Z  -d  sin  /3 
' min 

j2  {v  - v . )( e - e . ) 

max  min  max  min 


If  B < 77  / 2,  then  the  first  interval  is  (d  cos  B / Vmax  )<  t <_  min  [ (d  cos/?/  Vmjn  ) , 


(d  co  s cv/  V ) I ‘incl 
max  J 
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If  B-  7T-/ 2,  the  first  interval  is  of  zero  length,  and  the  formula  for  the  second 
interval  should  be  used  at  0. 

If  max  I (d  cos/3/V  . ),  (d  cos/3  /V  )]  < t < min  [ (d  cos  a/V  . ), 

(d  cos  cv/V  )1,  then  the  extremal  values  of  0 in  the  second  interval  are  on 
max  1 

arc-cosine  curves  and 
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If  min  I (d  cos  n / V ) , ( cl  cos  a/  V ) 1 < r<  max  | (cl  cos/?  / V . ),(d  cos/3/V  ) ] > 

lv  min  max  J-  ~ 1 r mm  max 

then  the  extremal  values  of  0 are  on  horizontal  section  boundaries  and 
. . cl  sin  /?  - sin  o 

P(T)  = ~zW  -v~Yo'  ~<r :>  - 

t max  min  max  min 

If  the  limits  for  the  second  interval  are  equal,  then  only  the  formulas  for  the 
first  and  third  intervals  should  be  used. 

If  nl  2<«  , then  the  third  interval  is  max[(d  cos/J/V  ), 

u max 

( d cos  (x/  V . ) 1 < r < ( d cos  a/V  ) and 
min  J ~ ~ max 
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If  <v<  7T /2,  then  the  third  interval  is 

max  I (d  cos/3/V  . ),  (dcosa/V  )]  < t < (d  cos  «/ V ) and 
L min  max  nun 
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If  rv  = Tj  / 2,  the  third  interval  is  of  zero  length,  and  the  formula  for  the  second 
interval  should  be  used  at  r = 0, 

Figures  V-2  and  V-3  present,  for  the  uniform  rectangular  velocity- 
azimuth  distribution,  the  time-lag  probability  density  function  p(r)  at 


U 

V /V  =0.1  and  V /V  = 0.9,  respectively,  for  various  30  sectors, 
min  max  min  max 

At  V /V  = 0.  1,  there  are  sharp  peaks  near  r = 0 in  Figure  V-2  (in 
min  max 

sharp  contrast  to  the  corresponding  inverse  velocity  model  of  Figure  III-6). 


At  V /V  = 0.9,  however,  there  is  a strong  resemblance  in  Figure  V-3 
min  max 


to  the  comparable  inverse  velocity  model  time-lag  probability  density  functions 
in  Figure  III  - 7 . 
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FIGURE  V - 2 


TIME-LAG  PROBABILITY  DENSITY  FUNCTION  FOR  THE  UNIFORM 

RECTANGULAR  VELOCITY -AZIMUTH  MODEL  AT  V /V  =0.1 

o min  max 

FOR  VARIOUS  30  SFCTOFS 
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FIGURE  V - 3 

TIME-LAG  PROBABILITY  DENSITY  FUNCTION  FOR  THE  UNIFORM 

RECTANGULAR  VELOCITY -A  ZIMUTH  MODEL  AT 

V . /V  =0.9  FOR  VARIOUS  30°  SECTORS 
min  max 


c. 


TAPERED-VELOCITY  UNIFORM -A  ZIMU  TH  DISTRIBUTION 


In  this  model,  the  probability  density  is  nonzero  within  a rectangle 
in  velocity -azimuth  space  and  is  the  product  of  the  velocity  probability  density 
function  p(V)  and  the  arrival  angle  probability  density  function  p(f))r 


p(V,  0)  p(  V)  p(  () ) . 

The  velocity  probability  density  function  p(V)  is  smooth  and  falls  off  to  zero 
at  the  minimum  and  maximum  velocities.  It  is  specified  over  three  equally 
long  intervals  by  quadratic  functions.  In  the  first  interval,  where 
V :LV<  (2V  . + V )/ 3, 

nnn  min  max 
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In  the  second  interval,  in  which  (2V  ^ + ^ max^  ^ ^ — ^niin  + 
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where  V = (V  + V )/2  is  the  mean  velocity, 
med  min  max 

where  (V  . +2V  )/3<LV<_V  , 

min  max  max 
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In  the  third  interval  , 
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Figure  V-4  shows  the  velocity  probability  density  function  p(V).  If  V and 

3. 

V designate  the  minimum  and  maximum  velocity,  respectively,  for  the  i-th 
b 

interval,  then 


V = 
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(3-i)  V . + i V 

mm  max 


and 


Vb  = 


(4-i)  V . + (i-1)  V 

min  max 


The  arrival -angle  probability  density  function  p(  9)  is  uniform  between 


0 . and  9 '■ 

mm  max 


- if  9 . <e<e 

tf  - f)  . min  max 


p(  9)  - \ max  min 


. i 


0 otherwise. 

As  in  earlier  models,  breaking  up  the  angular  range  Q . <9<9 

mm  max 

into  sections  in  accordance  with  the  procedure  depicted  in  Figure  III  - 4 and 
subsequently  adding  the  contributions  from  the  appropriate  sections  can 
simplify  determination  of  the  time-lag  probability  density  function  p(r). 

As  before,  « and  ft  denote  the  minimum  angle  and  the  maximum  angle,  re- 
spectively, of  each  section  after  mapping  or?  ■>  the  first  and  second  quadrant. 
At  the  time  lag  r , the  contribution  to  p(r)  from  each  section  is 

io 
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p(  V, 0 ) cos  9 d 9 , 


min 


where  <•>  and  a>  are  the  minimum  and  maximum  angles,  respectively, 
min  max 

at  the  time  lag  r after  the  rectangle  described  by  <y&.9£.ft  and  V . <LV <*V 

6 ° j r mm  max 

is  transformed  into  ( t,  0)  space.  The  term  |t|  in  the  denominator  may  be 
replaced  by  t if  the  integration  limits  are  reversed  when  t is  negative. 
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FIGURE  V-4 


TAPERED  VELOCITY  PROBABILITY  DENSITY  FUNCTION 


The  contribution  to  p(r)  is 
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in  the  first  velocity  interval, 
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in  the  third  velocity  interval,  where  ^ and  are  yet  to  be  specified. 

In  the  discussion  which  follows,  is  the  minimum  velocity  of  the 
i-th  velocity  interval  and  Vb  the  maximum  velocity.  If  n/Z<  a ^ li,  there  is 
a contribution  to  p(  r ) when  (d  cos  /*/Va)<r<(d  cos  ft  /V^,  In  this  case,  tin 
upper  integration  limit  is 
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when  (d  cos  /3/V  )<T<(d  cosn/V^. 

If  a<(i<nl  2,  there  is  a contribution  to  p(r)  when  (d  cos  /3/Vb)<r< 

(d  cosa/V  ).  In  this  case,  the  upper  integration  limit  is  dfa=  0 when 
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If  cn  <77-/  2 y3  and  r=  0,  there  is  no  simple  method  to  eliminate  the 

powers  of  r in  the  denominator  of  the  integrals  to  be  evaluated.  However, 

the  probability  P(  0 < r £.r)  that  r lies  between  0 and  r is 
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when  r is  small  enough,  and  the  probability  density  function  p(r)  is 
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since  the  integral  expression  is  simply  the  average  velocity  for  a velocity 

probability  density  function  symmetrical  about  (V  . + V )/2. 

1 min  max 


For  the  taper ed -velocity  uniform-azimuth  distribution,  Figures  V-5 

and  V-6  depict  the  time-lag  probability  density  function  p(r)  at  V . /V 

1 min  max 

0.  1 and  V /V  = 0.9,  respectively,  for  specified  30  sectors.  In  this 
min  max 

particular  model,  the  functions  p(")  are  smooth  and  continuous  unless 

I ) or  0 is  an  odd  integer  multiple  of  a/2.  The  probability  density 

min  max 

functions  p(i)  in  Figure  V-6  resemble  those  of  Figure  III-7,  IV-2,  and  V-3: 

when  the  minimum  and  maximum  velocities  arc  nearly  equal,  the  dominant 

influence  on  p(  r ) is  apparently  the  .pccific  angular  range  ()  . < 0 <0 

min  max 

of  the  uniform  azimuthal  distribution. 
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FIGURE  V-5 

TIME-LAG  PROBABILITY  DENSITY  FUNCTION  FOR  THE  TAPERED- 

VELOCITY  UNIFORM-AZIMUTH  MODEL  AT 

V /V  =0.1  FOR  VARIOUS  30°  SECTORS 
min  max 


} 


V-19 


P(t) 


() 

2 


P(T) 


min 


FIGURE  V-6 

TIME-LAG  PROBABILITY  DENSITY  FUNCTION  FOR  THE  TAPERED 
VELOCITY  UNIFORM-AZIMUTH  MODEL  AT 
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SUMMARY 

This  report  has  presented  two  algorithms  for  implementing  Wiener 
adaptive  multichannel  filtering  with  distributed  signal  models.  These  algo- 
rithms may  prove  useful  in  eliminating  or  ameliorating  problems  previously 
encountered  with  maximum  likelihood  adaptive  multichannel  filtering.  Among 
these  problems  are  mutual  cancellation  of  interfering  events,  signal  distor- 
tion, and  sensitivity  to  slight  deviations  from  an  ideal  plane-wave  signal 
model.  The  two  algorithms  discussed  require  some  method  for  estimating 
the  crosscorrelation  functions  between  the  signal  and  the  channels  entering 
the  adaptive  beamformer.  These  crosscorrelation  functions  are  estimated 
by  convolving  the  time-lag  probability  density  function  p(r)  corresponding  to 
a specified  velocity-azimuth  incoming-energy  distribution  with  the  signal 
autocorrelation  function,  which  is  approximated  by  averaging  the  input-channel 
autocorrelation  functions. 

Estimating  the  time-lag  probability  density  function  p(r)  for  various 
directionally-distributed  signal  models  is  an  interesting  problem  in  its  own 
right.  In  this  report,  three  basic  models  are  described.  The  first  is  an 
inverse  velocity  space  model,  in  which  the?  signal -model  probability  distri- 
butions arc  specified  as  a function  of  the  two-dimensional  inverse  velocity 
vector  u = V/(V  * V),  where  V is  the  incoming  energy's  apparent  velocity 
in  the  plane  of  the  array.  The  second  is  a distributed  ring  model  where  all 
incoming  energy  is  concentrated  at  a single  apparent  velocity  V with  respect 
to  the  plane  of  the  array,  but  is  distributed  over  a range  of  azimuths  according 
to  a known  probability  density  function  p^(0).  The  final  model  is  a velocity- 
azimuth  space  model,  in  which  the  signal  model  is  specific  d both  as  a function 
of  the  apparent  velocity  V in  the  plane  of  the  array  and  as  a function  of  the 
arrival  azimuth  0 . 
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